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£\j 1 We present predictions for multiplicities and single inclusive particle production in proton-lead 

y—{ , collisions at the LHC. The main dynamical input in our calculations is the use of solutions of the 

running coupling Balitsky-Kovchegov equation tested in e+p data. These are incorporated into a 
£N| realistic model for the nuclear geometry including fluctuations of the nucleon configurations. Particle 

— _ production is computed via either fct-factorization or the hybrid formalisms to obtain spectra and 

yields in the central and forward rapidity regions, respectively. These baseline predictions will be 

■ useful for testing our current understanding of the dynamics of very strong color fields against 
upcoming LHC data. 

m: 

! I I. INTRODUCTION 

9-; Heavy-ion collision programs performed previously at the RHIC and SPS accelerators have always benefited 

I tremendously from insight obtained by proton (deuteron)-nucleus collisions. It is expected, also, that the planned 
D . p+Pb run at the LHC will provide crucial benchmarks for the characterization of the hot QCD medium produced in 

■ lead-lead collisions, arguably a Quark Gluon Plasma (QGP). The absence of intricate final state effects induced by 
the formation of a QGP renders proton-nucleus collisions an excellent laboratory to study the so called initial state 

■ effects. These originate from the fact that the colliding nuclei do not behave as a mere incoherent superposition of 
1 their constituent nucleons. Rather, coherence effects are important and modify not only the partonic flux into the 

collision, but also the underlying dynamics of particle production in the scattering processes. A careful distinction 
CD of initial from final state effects is of vital importance for a proper characterization of the latter, as they may lead 

sometimes to qualitatively similar phenomena in observables of interest. Physics prospects for the p+Pb at the 
LHC run were recently compiled in 

Besides its role as a reference experiment, the p+Pb run at the LHC will provide access to kinematic regions 
never explored so far in nuclear collisions and thus carries great potential for discovery of new QCD phenomena 
on its own. In particular, the huge leap forward in collision energy with respect to previous high energy electron- 
nucleus or proton-nucleus experiments 1 will probe the nuclear wave function at values of Bjorken-x smaller than 
ever before. It is theoretically well established that at small enough values of Bjorken-x QCD enters a novel regime 
governed by large gluon densities and non- linear coherence phenomena @. The Color Glass Condensate (CGC) 
effective theory provides a consistent framework to study QCD scattering at small-x or high collision energies (for 
a review see e.g. 0,3)- It is based on three main physical ingredients: First, high gluon densities correspond 
to strong classical fields (f], Q , which permit ab-initio first principles calculation of "wave functions" at small x 
through classical techniques. Next, quantum corrections are incorporated via non-linear renormalization group 
equations such as the B-JIMWLK hierarchy or, in the large- N c limit, the BK equation 0, [|| which describe the 
evolution of the hadron wave function towards small x. The non-linear, density-dependent terms in the CGC 
evolution equations are ultimately related to unitarity of the theory and, in the appropriate frame and gauge, 
can be interpreted as due to gluon recombination processes that tame or saturate the growth of gluon densities 
for modes with transverse momenta below a dynamically generated scale known as the saturation scale, Q s (x). 
Finally, the presence of strong color fields A ~ 1/g leads to breakdown of standard perturbative techniques to 
describe particle production processes based on a series expansion in powers of the strong coupling g. Terms of 
order gA ~ C(l) need to be resummcd to all orders. The CGC provides the tools to perform such resummation 
although the precise prescription for the resummation may vary from process to process or colliding system. 

An important step towards promoting the CGC to a practical phenomcnological tool has been performed 
recently through the calculation of higher order corrections to the BK-JIMWKL equations [9HL^|. In particular, 
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1 The expected center of mass collision energy for the p+Pb run is 5 TeV, to be compared to a maximal energy of 200 GeV for RHIC 
d+Au collisions. 
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the calculation of running coupling corrections to the evolution kernel of the BK equation (referred to as rcBK 
henceforth) made it possible to describe various data at high energies in terms of solutions of the rcBK equation, 
thus closing the gap between first principles theory calculations and data fl3j . 

Despite the fact that CGC effects arc expected to be enhanced in nuclei versus protons, which is due to the larger 
valence charge densities per unit transverse area in nuclei, so far the most exhaustive searches of the saturation 
phenomenon have been performed using data on proton reactions. This is mainly due to the large body of high 
quality DIS data on protons at small-x. Thus, the global fits performed by the AAMQS collaboration in |14l-ll6| 
demostrated that the rcBK equation successfully accounts for the x- dependence of inclusive structure functions 
measured in electron-proton collisions at HERA at small- a; (see also [17| for fits to diffractivc data). Here, we 
shall use the resulting AAMQS parametrizations of the dipole-proton scattering amplitude as a basic building 
block for the nuclear unintcgratcd distributions. Furthermore, as discussed in more detail below, single inclusive 
p± spectra in proton+proton collisions at Tevatron, LHC and forward RHIC data are well described in the CGC 
framework [L8l - r2C| . In case of heavy- ion collisions, the most compelling indications of CGC effects in data are 
due to suppression phenomena observed at RHIC d+Au collisions in the forward rapidity region. Actually, single 
particle distributions in both p+p and d+Au collisions and the depletion of nuclear modification factors RdAu 
were well described by CGC calculations in [IHOjl. However, forward RHIC data are close to the kinematic limit 
of phase space, probing the projectile wave function at large values of x — > 1 and the fragmentation functions at 
z — > 1. This allows for alternative descriptions of the data based on large- a; phenomena such as energy loss [2l[ . A 
clearer signal of gluon saturation is provided by the observed disappearance of back-to-back di-hadron correlations 
in d+Au collisions. CGC calculations predicted first (22[, and later on quantitatively accounted for [23|, [24[, such 
suppression. The disappearance of the away-side peak can be understood in terms of multiple scatterings controlled 
by a dynamical transverse scale, the saturation scale of the nucleus. Even though no description of such data that 
does not invoke saturation effects is known to date (see however [HI]), it has been argued that large- a; effects like 
energy loss or enhancement of double scatte ring processes may blur the interpretation of angular decorrelation 
between hadron pairs as due to CGC effects [26[ . Additionally, the good description of the energy and centrality 
dependence of integrated multiplicities in RHIC Au+Au and d+Au collisions and in Pb+Pb collisions at the LHC 
in CGC models lends further support to the idea that saturation effects are relevant in present data, although the 
larger degree of modeling involved in their calculation prevents making any definite claims. In summary, although 
it may be argued that the present data do not provide sharp and direct evidence for neither the saturation regime 
of QCD nor the CGC approach to saturation, it is fair to say that HERA, RHIC and available LHC data lead 
to a rather coherent picture, with several phenomena finding their natural interpretation in terms of high density 
gluonic systems together with a consistent quantitative description in the CGC framework at its present degree of 
accuracy. 

Regardless of the question whether the CGC is the most suited framework for their description there is broad 
consensus that coherence effects are important for the interpretation of present data on heavy ion collisions. In fact 
most -if not all- of the different phenomenological approaches for the description of particle production -both in 
the soft or hard sector- include physical ingredients related to either the depiction of gluon content of nuclei (e.g. 
the leading twist shadowing of the nuclear parton distributions in the collincar factorization formalism or string 
fusion processes in Dual Parton or percolation models), the breakdown of independent particle production from 
the individual participant nucleons (e.g. through the presence of energy dependent cut-offs in event generators like 
HIJING or HYDJET) or effective resummations of multiple scatterings 2 . All these ingredients are akin, at least 
at a conceptual level, to those dynamically built in the CGC, although they are formulated in very different ways. 

The p+Pb run at the LHC will provide an excellent -and probably in the near future unique- possibility 
to disentangle the presently inconclusive situation on the role of CGC effects and also to distinguish among 
different approaches to describe high energy scattering in nuclear collisions. On the one hand, the LHC shall 
bring us closer to the limit of asymptotically high energy in which the CGC formalism is developed, thus reducing 
theoretical uncertainties on its applicability. Equivalcntly, the value of the saturation scale is expected to be a 
factor ~ 2 ~ 4 times larger than at RHIC, so saturation effects should be visible in a larger range of transverse 
momenta, deeper into the perturbative domain. On the other hand, the much extended energy reach of the LHC 
will allow measurements far from the kinematic limit up to very forward rapidities, thus minimizing the role of 
large- x effects which obscured the interpretation of forward RHIC data. 

Although more exclusive obscrvablcs like di-hadron or hadron-photon correlations are expected to better dis- 
criminate between different approaches, a first test for models of particle production shall come from data on 



2 Here we just mention a few well known examples; a rather exhaustive compilation of phenomenological works for the description of 
particle production in HIC can be found in e.g 12711 . 
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inclusive multiplicities and single particle distributions as they are much easier to obtain experimentally. In this 
work we aim to provide up-to-date predictions for these observables within the framework of the running coupling 
BK Monte Carlo (rcBK-MC) previously presented in [28| and built as an upgrade of the original KLN-MC code 
[2^]. The rcBK-MC model presented here attempts to incorporate up-to-date tools on the CGC theoretical side, 
namely the use of solutions of the rcBK equation as main dynamical input to describe the x and transverse momen- 
tum dependence of nuclear unintegrated gluon distributions, combined with the empiric information gained from 
the analysis of e+p, p+p and RHIC nuclear collisions. The Monte Carlo treatment allows for a realistic treatment 
of the collision geometry including in particular realistic nuclear density distributions as well as fluctuations of 
the configurations of nucleons. This is important in order to obtain the number of collisions N co u used in the 
definition of the nuclear modification factor R p a consistently with the computed p+A spectra, within one single 
framework. Also, fluctuations of the number of target nucleons are particularly relevant for particle production in 
the high density regime where the number of small- x gluons is not linearly proportional to the number of valence 
charges. 

There are two distinct but related approaches to hadron production in high energy asymmetric (such as proton- 
nucleus or very forward proton-proton and nucleus-nucleus) collisions. In such collisions, particle production 
processes in the central rapidity region probe the wave functions of both projectile and target at small values of 
x. Here, one may employ the fc t -factorization formalism where both the projectile and target are characterized in 
terms of their rcBK evolved unintegrated gluon distributions (UGDs). 

At more forward rapidities, on the other hand, the proton is probed at larger values of x while the target nucleus 
is shifted deeper into the small- a; regime. Here, fct-factorization fails to grasp the dominant contribution to the 
scattering process. Rather, the hybrid formalism proposed in ref. [30[ and embedded recently in a Monte Carlo 
treatment of the nuclear geometry (l9l l20j is more appropriate. In the hybrid formalism the large- X degrees of 
freedom of the proton are described in terms of usual parton distribution functions (PDFs) of collinear factorization 
which satisfy the momentum sum rule exactly and which exhibit a scale dependence given by the DGLAP evolution 
equations. On the other hand, the small-x glue of the nucleus is still described in terms of its UGD. Recently the 
hybrid formalism has been improved through the calculation of inelastic contributions that may become important 
at high transverse momentum [3l|. We shall investigate the effect of this new production channel in our predictions. 
Unfortunately, no smooth matching between the k t -factorization and hybrid formalisms is known to date. Also, 
their corresponding limits of applicability -equivalently the precise value of x at which one should switch from one 
to the other- have only been estimated on an empirical basis. We shall explore the stability of CGC predictions 
using one or the other formalism as one varies the kinematics of the detected hadron. 

One important goal of this work is to systematically assess the model (and implementation) uncertainties. Such 
calibration is mandatory in order to actually prove or disprove different approaches by comparison to data. This 
is specially so in the case of semi-inclusive observables such as the ones discussed in this work where "only" 
quantitative differences arise between predictions from different approaches (contrary to the case of more exclusive 
observables). 

In the rcBK-MC model presented here, uncertainties originate from either the modeling of non-pcrturbative as- 
pects of the collision or from high-order corrections to fixed order perturbative calculations for particle production. 
Among the former are the initial conditions for rcBK evolution, which are only partially constrained by global fits 
to c+p data, or modeling of the impact parameter dependence of the nuclear UGDs. Perturbative tools such as 
the rcBK equation are ill-suited to provide a simultaneous, unified description of both the Bjorken-x and impact 
parameter dependence of the UGD of a nucleus (or any other hadron) since the latter is related to the physics of 
confinement. This makes some degree of modeling unavoidable. Concerning the uncertainties due to fixed order 
perturbative calculations, they are either absorbed in if-factors or explored through the variation of factorization 
and running coupling scales. The if -factors discussed here should also absorb effects not included in the CGC 
approach. 

II. THE NUCLEAR UNINTEGRATED GLUON DISTRIBUTIONS 

The paucity of small- a; data with nuclei prevents a direct empiric determination of the nuclear unintegrated 
gluon distribution at small-x. Here we shall build the nuclear UGD entirely from the nucleon UGD extracted by 
the AAMQS collaboration through global fits to inclusive c+p data on structure functions. 

Before proceeding further, let us briefly recall the basic features of the AAMQS approach. The main dynamical 
ingredient is the dipole formulation of deep inelastic scattering (DIS). The AAMQS fits rely upon is the qq dipolc 
scattering amplitude off a hadron, A/f(»", x, R), where x is the usual Bjorken scaling variable in DIS, r is the dipolc 
transverse size and R the transverse position at which the hadron target is probed. The index F refers to the fact 
that the dipole constituents -quark and antiquark- belong to the fundamental representation of SU(3). In the 
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large i\T c -limit the B-JIMWLK equations reduce to a single, closed equation for the x-depcndcnce of the dipolc 
amplitude, which is the so-called BK equation. Under the assumption of translational invariancc implicit in the 
AAMQS approach, and hence omitting the R-dcpendcncc of the dipolc amplitude accordingly, the BK equation 
equation including running coupling corrections (referred to as rcBK in what follows) reads 



dM F (r,x) 



d \ti{xq/x) 
where r = r 1 + r 2 3 and K7 



d 2 Ll K run (r,r 1 ,r 2 ) [A/>(n,x) + N F (r 2 ,x) - M F {r,x) - M F {r 1 ,x)M F {r 2 ,x)] 



is the evolution kernel including running coupling corrections: 
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In practical implementations, the running coupling in Eq. ^ is regularized in the infrared by freezing it to a 
constant value a& = 0.7. 

Solving the BK equation is an initial value problem, i.e. it is well defined only after initial conditions at the 
initial evolution scale, Xq = 10~ 2 in the AAMQS fits, and for all values of the dipole size r have been provided. 
This introduces free parameters, ultimately of non-perturbative origin, to be fitted to data. In the AAMQS rcBK 
fits to HERA data the initial conditions are taken in the form 



N F (r,x=xa) = 1 — cxp 



(r 2 Q 2 
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where A = 0.241 GeV, Q 2 proton is the saturation scale at the initial scale xq and 7 is a dimensionless parameter 
that controls the steepness of the tail of the unintegrated gluon distribution for momenta above the saturation scale 
h > Qso, proton- Both Q^o, proton and 7 are fitted to data. Although the AAMQS fits clearly favor values 7 > 1, 
they do not uniquely determine its optimal value (and neither does the analysis of forward RHIC data performed 
in ref. [l9j|). Rather, different pairs of (Qgo.proton' 7) parameters provide comparable values of x 2 /d.o.f ~ 1. The 
reason for this behavior is that they are correlated with other parameters, such as the overall normalization of 
the 7*-p cross section, and also that HERA data is too inclusive to constrain exclusive features of the proton 
UGD. In order to account for such uncertainty we shall consider two of the AAMQS sets, corresponding to 
(<9?o, P roton:7) = (0-168 GeV 2 1.119) and (0.157 GeV 2 , 1.101). Additionally, we shall also consider the McLerran 
Venugopalan (MV) model [H Q which corresponds to Eq. ([3]) with 7 = 1; the MV model is well established 
theoretically (for very large nuclei, gA 1 ^ 3 3> 1) and has been used frequently in the literature. We note that 7 > 1 
for the proton may arise from corrections to the effective action of the MV model of higher order in the valence 
color charge density [32I ] . Such corrections are expected to decrease with increasing nuclear thickness. Therefore, 
it is conceivable that the dipole-nucleus scattering amplitude may be better represented by the MV model than by 
initial conditions with 7 > 1. However, using the AAMQS (7 > 1) i.e. for the proton simultaneously with an MV 
i.e. for the nucleus leads to a monotonical increase (with p t ) of the nuclear modification factor i? p p D which is due 
to the different power-law fall off of the respective UGDs. We shall not study this option in detail in this work. 
Also, we recall that in the AAMQS fits the strong coupling is evaluated according to the following expression: 
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where j3 = 11 — §iV/, Nf = 3, A = 0.241 GeV and the constant C under the logarithm accounts for the 
uncertainty inherent to the Fourier transform from momentum space, where the original calculation was performed. 
A parameter \i is introduced to regulate the strong coupling for large dipole sizes, and fixed by the condition 
a s {po) = a/ r . The (<5 2 proton , 7, a/ r , C)-valucs considered here are shown in TablcQ] 



UGD Set 


QsO, proton (GeV 2 ) 


7 


a fr 


G 


MV 


0.2 


1 


0.5 


1 


gl.119 


0.168 


1.119 


1.0 


2.47 


gl.101 


0.157 


1.101 


0.8 


1 



TABLE I: Summary of parameters for the three dipole-proton scattering amplitudes (or UGDs) considered in this work. 



3 We use the notation v = |v| for two-dimensional vectors throughout the paper. 
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In the large- N c limit the gluon dipole scattering amplitude required for the unintegrated gluon distributions can 
be obtained from the quark dipole scattering amplitude that solves the rcBK equation: 

Na{t, x, R) = 2A/>(r,.T,R) - A/|(r,x,R) (5) 

where the subscript A simply refers to the fact that gluons belong to the adjoint representation of SU(3). Note 
that this relation entails that the saturation momentum relevant for gluon scattering is larger than that for quark 
scattering by about a factor of 2, at the initial rapidity x — x$. 

The nuclear UGDs needed for particle production in the fc t -factorization and hybrid frameworks are related to 
the quark and gluon dipole scattering amplitudes via 2-dimensional Fourier transforms. In particular, the UGD 
entering the fc t -factorization formula is given by 

<p{k,x,R) = ^ f d 2 ver^V 2 r Af A (r,x,R). (6) 
a s {k) (2^ J 

The function tp is dimcnsionless and corresponds to the number of gluons per unit transverse area and per transverse 
momentum space cell. 

Eq. ^ was written originally for fixed coupling. In order to be consistent with our treatment of the small- a; 
evolution, we have extended it by allowing the coupling in the denominator to run with the momentum scale. 
As explained in (28|, this modification turns out to be important to obtain a good description of the centrality 
dependence of the charged hadron multiplicities measured in Au+Au and Pb+Pb collisions at RHIC and the 
LHC respectively, which are otherwise too flat. In turn, the following UGDs in the fundamental and adjoint 
representations are needed in the hybrid framework: 

%) (*, x,TL) = Jd 2 r e" lkr [l - M F[A) (r, x, R)] . (7) 

Note that although the AAMQS approach assumes translational invariance -i.e. impact parameter independence— 
of the dipole scattering amplitude over transverse distances of the order of the nucleon radius, i?/v, we have made 
explicit such dependence in equations Eqs. ([5][7]) for consistency with the notation employed for nuclear UGDs 
where the impact parameter dependence must be considered. 

Let us now discuss our model for the nuclear UGD, starting from the one for a nucleon 4 . First, we assume 
that the functional form of the quark dipole scattering amplitude off a nucleus at the initial saturation scale 
x = xq = 0.01 is the same as for the quark dipole scattering amplitude off a proton but with a shifted initial 
saturation scale that depends on the local density at every point in the transverse plane, R. In other words, we 
shall replace 

QsO, proton ~~ ^ QsO, nucleus (R) (8) 

in Eq. ([3]) in order to define the initial conditions for the evolution of the quark dipole scattering off a nucleus at 
every transverse point in the nucleus. Then, the initial conditions given by Eqs. ((3j) and ((8j) arc evolved locally using 
the impact parameter independent rcBK evolution defined by Eqs. ([TJ2]). This provides the full (r, x)-dependence 
of the quark dipole-nucleus scattering amplitude at every point in the transverse plane of the nucleus. Finally, 
Eqs. ([HE]) are used to calculate their Fourier transform which provide the complete transverse momentum fc t , 
Bjorken-a; and impact parameter (R) dependence of the nuclear UGDs. 

To complete our discussion of the initial conditions we explain how we construct Q 2 nuc i ous (R)- We treat the 
transverse positions of nucleons as random variables following a two-dimensional projection of the Woods-Saxon 

distribution, Ta(R). Each configuration consist of a list of random coordinates r,, i = 1 A, for the locations of 

the different nucleons in the transverse plane; A denotes the atomic mass number of the nucleus. Multi-nucleon 
correlations are neglected except for imposing a short-distance hard core repulsion which enforces a minimal 
distance rj 0.4 fm between any two nucleons. 

Every such configuration defines a different local density in the transverse plane of the nucleus. Obviously, the 
smallest non-zero local density corresponds to the presence of a single nucleon. On the other hand, rare fluctuations 
where a large number of nucleons is encountered at a given transverse position can occur. Such configurations 
correspond to a high initial saturation scale, Q 2 Q nucleus (R). For a given configuration, the initial saturation 



4 Strictly speaking, the AAMQS fits provide information only on the proton UGD or dipole scattering amplitude. We shall assume 
that isospin effects are negligible for the gluon distributions at small- ie and consider it equal to the one of a neutron. 
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momentum Q 2 s0 nuc i 0US (R) at the transverse coordinate R is taken to be proportional to the local variance of the 
density of valence charges which is itself proportional to the number of overlapping nuclcons at that transverse 
point, A(R): 

a ) QsO,nuclcus(R) = N ( R ) Qlo, proton • ( 9 ) 

We shall refer to this option as natural in what follows. Note that these valence charge sources generate the 
small-x gluon fields described by the UGD coherently. For the AAMQS initial condition, the prescription ([9j may 
lead to an inconsistent definition of the nuclear modification factor in the limit of large transverse momentum 
because additivity of the dipole scattering amplitude in the number of nucleons, Af(r; xq, R) ~ Ta(R) at small r, 
is violated. To fix this issue we shall also consider a second possibility to generate Q 2 s0 nuc i eus (R): 

b) Q?0,nucleus(R) = N(R) lh Q 2 s0< proton , (10) 

where 7 is one of the parameters in the initial condition, see Eq. This ansatz -to which we shall refer as 
modified— is motivated by the requirement that nuclear modification factors should go to unity at large transverse 
momentum, i? p +p D (pt 3> Q s ) — > 1, as we shall discuss in detail later. While we presently lack a solid theoretical 
derivation for Eq. (flUf we view it as a phenomenological way to ensure that the nuclear UGD is additive in the 
number of nucleons at small dipole sizes r or high intrinsic transverse momenta k t . 

In both cases, the total number of nucleons that overlap at a given transverse point, 7V(R), is calculated via a 
simple geometric critcrium: 

7V(R) = ^e(J^-|R-r i |) . (11) 

Some care must be exercised in choosing the transverse area oq of the large- a; partons of a nucleon. Q s q, proton 
corresponds to the density of large- a; sources with x > xq and should therefore be energy independent (recoil of 
the sources is neglected in the small- a; approximation). We therefore take ao — 42 mb to be given by the inelastic 
cross-section at = 200 GcV. However, cto should not be confused with the energy dependent inelastic cross 
section <7i n (s) of a nucleon which grows due to the emission of small- a; gluons. 




1 10 1 10 

k, (GeV/c) k, (GeV/c) 



FIG. 1: Left: unintegrated gluon distributions tp(kt; x) for different values of the initial saturation scale evolved to a; = 3-10 
for MV initial conditions. Right: Fourier transform of the dipole scattering amplitude J\fA(k;x) for a single nucleon and 
two different initial conditions: MV 7=1 ' 119 (solid) and MV (dashed) at rapidities Y = ln(xo/x) = 1.5 and 6. 

In Fig. [T] we plot the UGD <p{kt\ x) for three different initial saturation scales at x = 3 ■ 10 -4 versus transverse 
momentum. The UGD corresponding to a single nucleon peaks at about kt — 1 GeV. The UGDs for larger Q 2 s0 
illustrate the shift corresponding to a 6-nucleon and 12-nucleon target, respectively. These curves illustrate why 
one hopes that particle production at high energies and/or for heavy targets would become insensitive to physics 
at the confinement scale while instead being dominated by the semi-hard scale generated dynamically. 

The plot on the right shows the different slopes for the two initial conditions (MV 7=1119 versus MV). The 
AAMQS UGD exhibits a substantial suppression of the high-fc t tail as compared to MV model initial conditions. 
At fixed kt the suppression diminishes with evolution to higher rapidities. We shall show below that the AAMQS 
UGD provides a much better description of semi-hard p t spectra in p+p collisions at high energies. 
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FIG. 2: Left: xi, 3 coverage at y = 2.2, 3.2 and 4.0 at RHIC energy s/s = 200 GeV. Hadron p±_ is limited to p T < 7.5 (3.5) 
GeV/c at y = 3.2 (4.0) at RHIC. Right: xi t 2 coverage at y — 4, 6 at LHC energy. 



III. PARTICLE PRODUCTION: fc t -FACTORIZATION AND HYBRID FORMALISMS 

In this section we summarize how particle production is calculated from the UGDs described above. It is useful 
to first recall some elementary kinematics. For inclusive production of a single parton with transverse momentum 
Pt and rapidity y, without detection of the recoiling particle(s) in the opposite hemisphere, the 2 — > 1 kinematics 
is such that projectile and target fields are probed at light cone momenta = (Pt/\/s)e ±v . For the upcoming 
p+Pb run at the LHC we shall assume a collision energy of y/s = 5 TeV. Hence, the production of particles with 
transverse momentum p t < 20 GeV in the central rapidity region \y\ < 1 probe both the projectile and target wave 
functions at small values of x below our initial condition at Xq = 0.01 5 . In this case we shall use the so-called 
fc t -factorization approach. On the other hand, at more forward rapidities y > 2, towards the proton fragmentation 
region, the momentum carried by the projectile parton grows large and so we shall resort to the hybrid formalism. 

To outline the limits of applicability of our model let us first mention that in our terminology "small- a;" refers to 
x < xq = 0.01. This is the limit of applicability of the AAMQS fits to HERA data and coincides with parametric 
estimates for the validity of the CGC approach and of coherent interactions. To estimate the range of pt where 
our calculations might be valid we again resort to the AAMQS fits: they start exhibiting some tension with the 
data when extended beyond Q ~ 7 -f- 10 GeV. This may indicate failure of the CGC approach based on rcBK 
rcsummation for higher virtualities or transverse momenta. These limits of applicability are, of course, merely 
indicative, since one does not expect a sharp boundary but rather a smooth transition. The kinematic window 
accessible at RHIC and LHC, respectively, is illustrated in Fig. [5J 



A. fct-factorization 

According to the fc t -factorization formalism [33j , the number of gluons produced per unit rapidity at a transverse 
position R in A+B collisions is given by 

dN A+B^g ^ 1 d(T A+B^g 

dy cPpt <PR ~V s dy d 2 Pt d 2 R ' (1) 
where <j s represents the effective interaction area and a A+B ^ 9 is the cross section for inclusive gluon production: 



dyd? Pt d 2 R C F p\ 



Kk T^^I V / d 2 ba s (Q)<p P l^^, Xl ;b U [^L—^^ 2 -R^b . (13) 



5 A more precise estimate should consider the kinematic shift induced by the convolution of the primary produced parton with the 
fragmentation function, which is indeed taken into account in our calculation later on. 
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y and pt are the rapidity and transverse momentum of the produced gluon, respectively, while x\ t i = 
(Pt/\/sNN) exp(±y) and Cf = (N 2 — 1)/2N C . As noted before, wc assume that the local density in each nucleus 
is homogeneous over transverse distances of the order of the nucleon radius Rn- Thus, the 6-integral in Eq. (|13l) 
yields a geometric factor proportional to the transverse "area" of a nucleon which cancels with a similar factor 
implicit in a s from Eq. (|12[) . modulo subtleties in the definition of <r s . Note, also, that (fTU| is symmetric under 
projectile <-> target exchange if, simultaneously, one lets y — > —y. 

Eq. (|13p was written originally for fixed coupling. For consistency with our running coupling treatment of the 
small- a; evolution of the UGDs, we allow the coupling to run with the momentum scale. The argument of the 
running coupling in Eq. (fl3]) is chosen to be Q — max{|p t + fc t |/2, \p t — k t \/2}. However, we found rather weak 
sensitivity to the particular choice of scale because ip — > as k t — > due to the saturation of Af(r) at large 
dipolc sizes r, see above. In principle one could improve on this educated guess by using the "running coupling 
fc t -factorization" formula derived recently in ref. [H| . Most importantly, the ^-dependence of the dipole scattering 
amplitude obtained by solving the rcBK equation encodes all the collision energy and rapidity dependence of the 
gluon production formula Eq. (|13[) . 

The normalization factor K k introduced in the fct-factorization formula (|13|) above lumps together higher-order 
corrections, sea-quark contributions and, effectively, other dynamical effects not included in the CGC formulation. 
It can be fixed approximately from the charged particle transverse momentum distribution in p+p collisions at 
7 TeV, see below. Although its precise value depends on the UGD and on the fragmentation function, we typically 
find K k ~ 1.5-3, which appears reasonable. 

Eq. (flU|) is the starting point for all observables discussed below. In particular, the charged particle multiplicity 
and the transverse energy can be obtained by integrating over the transverse plane and p t , 

dN ch 2 f , 2 „ f - dN A+B ^ 



dy . ^JfRjfpt——, (14) 
dE t f f , dN A+B ^s 



Note that a low-p t cutoff is not required since the integration over k t in (|13|1 extends only up to p t . The saturation 
of the gluon distribution functions guarantees that the dominant scale in the transverse momentum integrations 
is the saturation momentum. Similar to other CGC-based approaches, our Eq. (|14p assumes that the total hadron 
multiplicity is proportional to the initial gluon multiplicity through an (energy and centrality independent) gluon 
multiplication factor n g 6 . In order to reproduce both RHIC and LHC data on charged hadron multiplicities in 
heavy-ion collisions we fix it to be n g ~ 5; small adjustments of this normalization factor may be required for p+p 
and p+Pb collisions as discussed below. This could be due to the fact that we here assume local rcBK evolution 
without explicit impact parameter dependence; further, due to our ignorance about how to hadronize the small- a: 
gluons (which also determines the y — > rj Jacobian given below) etc. The precise value of n g does of course depend 
on the value of the if -factor which has been fixed independently since dN/dy only involves their product. 

The upper limit in the integrals over the gluon transverse momentum in dN^/dy and dE t /dy has been taken as 
pmux _ yi GeV; if the integrals are extended further then a slight adjustment of the normalization factors n g and 
K may be required. In order to compare our results for initial gluon production to the final state distributions 
of charged particles one has to translate the rapidity distributions into pseudo-rapidity distributions through the 
y — >• 1] Jacobian 

dN c h cosh?? dN c h dE t coshf? dE t 

d, l ~ Jcosh 2 n + myp2 d y ' dr > ~ J 'cosh 2 T, + m?/P* dy ' 



with y = iln(ycosh ri + m 2 /P' 2 + sinh?y)/(y cosh T] + m 2 /P 2 — sinhr/). Wc assume that in this Jacobian 

m = 350 McV and P = 0.13 GeV + 0.32 GeV(Vs/l TcV) 0115 which leads to a reasonably good description of 
the pseudo-rapidity distribution of charged particles in p+p collisions at LHC energies, see below. 

In turn, the single inclusive spectra at perturbatively large transverse momenta can be obtained by folding 
Eq. (|13[) for gluon production with the corresponding gluon fragmentation function: 

djyA+B-yhX r r j / \ jprA+B^g 



dyd 2 p t J J z 2 9 V kt ' J dyd 2 q t d 2 R 



This factor does not enter Eq. JTSJ for the transverse energy because we assume that final-state gluon showering and hadronization 
conserves the energy per unit rapidity. 
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In (|17[) the integral over the hadron momentum fraction is restricted to z > 0.05 to avoid a violation of the mo- 
mentum sum rule. The scale dependence of the FF of course emerges from a resummation of collincar singularities 
via the DGLAP equations and so its use in the fcj_-factorization formula is not entirely justified. 



B. Hybrid formalism 



Moving away from central rapidity towards the projectile fragmentation region its wave function is probed at 
larger and larger momentum fraction x\ which will eventually exceed x = 0.01. In this case the so-called hybrid 
formalism (30l | is better suited for particle production. We shall employ the following expression for the differential 
cross section for production of a hadron with transverse momentum k and pseudorapidity 7 r\: 



dNP 



A-^hX 



di] d 2 k 



= K' 





' dN h - 


+ 


' dN h ' 




( 


drj d 2 k 


el 


drj d 2 k 


incl^ 



(18) 



where the subscripts el and inel stand for elastic and inelastic contributions 8 , respectively. We again allow for the 
presence of a K- factor, K h , to absorb higher order corrections. The first term in Eq. (fT8|) . the clastic contribution, 
is given by [3(| 



[ dN h ' 


1 f 1 dz 


drj d 2 k 


e] - (270^, z 2 



L q 



+ x 1 f g/p (x 1 ,Q 2 ) N A (^2 A) D h/g (z,Q 2 ) 



(19) 



and corresponds to scattering of collinear partons from the projectile on the target. The 2—^1 kinematics sets 
x i,2 = (pt/ z^/snn) cxp(±y) and xf — (j>t/ \Jsnn) cxp rj. The projectile is described by standard collinear parton 
distribution functions (PDFs) but its partons acquire a large transverse momentum k due to (multiple) scattering 
from the small- a; fields of the nucleus which arc described by the corresponding UGDs in the adjoint or fundamental 
representation N^ F ^, see Eqs. 0. The hadronization of the scattered parton into a hadron is described by the 
usual fragmentation function (FF) of collinear factorization, D^/j- Both the PDF and the FF are evaluated at the 
factorization scale Q. We shall explore the sensitivity to the choice of factorization scale by letting it vary within 
the range Q = (k/2, 2k). 

The inelastic term in Eq. (|T8|) has been calculated recently in ref. [U . It reads 



dN h 



dr\ d 2 k 



a s (Q) 
2tt 2 



dz z 



d 2 q 
(2tt)' 



q 2 N F (x 2 , q)X] 



■ £ 

*,j=q,q,g 



"■, jar, ao/a 



)D h ,Az,Q 2 ),m 



where P^/j are the LO DGLAP splitting functions for the different parton species i, j — q, q,g. Note that endpoint 
singularities for q —> q and g — s- g splitting are regulated via the usual "+ prescription" ; therefore, the contribution 
from Eq. (|20p is actually negative in parts of phase space. Explicit expressions for the weight functions Wi/j(^) 
are given in Eqs. (74-77) of ref [3l[ and shall not be repeated here. 

The inelastic term corresponds to an alternative channel for hard production: partons with high transverse 
momentum can occur in the wave function of the incoming proton due to large-angle radiation. Those may then 
scatter off the target with only a small momentum transfer to finally fragment into a high-p t hadron. It should be 
noted that the inelastic term accounts for part of the full NLO corrections to the hybrid formalism. A calculation 
of the full NLO corrections to the hybrid formalism has been recently presented in [3^, [3(| . While a numerical 
implementation of the full NLO corrections would be necessary, such task is beyond the scope of this paper and 
we leave it for future work. The evaluation of the inelastic term in this work provides an estimate of the numerical 
importance of the full NLO corrections. 

The inelastic contribution involves an additional power of the coupling a s but also comes with a factor ~ 
\og(k 2 /Q 2 T ) [3lj and so is expected to be significant at high transverse momenta and not too forward rapidities, 
far from the kinematic boundary. These expectations have been first verified numerically in ref. j37j . Finally, we 
note that the scale for the running of the coupling in Eq. (f20)) cannot be determined from the calculation of ref. [3l| 



7 At forward rapidities the distinction between r\ and y becomes less relevant. 

8 As a NLO contribution, the latter need not be positive definite, see below. 
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-a full NNLO calculation would be necessary to determine the scale for the running of the coupling at NLO- and 
should be considered as a free parameter. In order to asses the uncertainty related to the choice of scale for the 
coupling in the inelastic term, we shall consider two possibilities: assuming a constant value a s =0.1 (similar to 
what was done in an earlier study f37| ) or, alternatively, assuming one- loop running at the factorization scale Q 
by replacing a s — > a s (Q) in Eq. (p?0"|) . 

IV. THE BASELINE: PROTON+PROTON COLLISIONS 

In this section we first present our results for proton-proton collisions. We restrict to CM energies in the TeV 
regime so that the typical parton momentum fractions involved in semi-hard hadron production remain small. 
A good description of p± distributions in p+p is of course required for a reliable calculation of the spectra in 
minimum-bias p+A collisions as well as for the R p a nuclear modification factor. Furthermore, it is important to 
check consistency of the UGD in DIS and hadronic collisions. In what follows we shall use three different sets of 
fragmentation functions: LO-KKP [H] and DSS [H,|40[ at L0 and NL0 > respectively. 




P t (GeV) p, (GeV) 



FIG. 3: Transverse momentum distribution of charged particles in the central region of P+P collisions at ^/s = 1.96 TeV 
(left), and p+p collisions at ^/ s — 7 TeV (right). CDF and CMS data from refs. [4l| and [43]; respectively. 

In Fig. [3]we show the transverse momentum distributions of charged particles in the regime of semi-hard pt for 
inclusive p+p collisions at y/s = 1.96 TeV and 7 TeV, respectively. We have tested various combinations of UGDs 
and FFs and in each case matched the A-factor to the data at p± = 1 GeV. 

First, we note that the UGD sets with the steeper initial gluon spectrum (UGD sets glll9 and gllOl) lead to 
significantly better agreement with the data than the classical "MV model" initial condition with 7 = 1; the latter 
leads to a p± spectrum far outside the experimental error bars. This establishes consistency within our framework 
of the UGDs with DIS and hadron-hadron collisions, since those two sets provide a much better x 2 /d.o.f in fits to 
e+p data than the MV one. Recall, also, that in hadronic collisions x ~ p± for spectra at fixed rapidity. Hence, 
the shape of the spectra does provide a direct test of the rcBK evolution speed. 

In Fig. [4] we check that changing the scale in the FF from Q = p±/2 to Q = p± to Q = 2p± is mainly absorbed 
into a redefinition of the A"- factor; in what follows we shall fix Q = p±. Similarly, switching from DSS-LO to 
DSS-NLO FFs essentially leads to identical spectra once the A-factor is adjusted; see Fig. [3] (left). On the other 
hand, we obtain slightly harder spectra with DSS versus KKP fragmentation functions, as expected. 

We now turn to the pseudo-rapidity dependence of pj_-integrated multiplicities. As already mentioned above, 
here we can not convolute the gluon spectrum with a fragmentation function. Instead, our crude "hadronization 
model" consist in assuming that on average over events the number of charged particles is proportional to the 
number of initially produced gluons. The number n g of final hadrons per initial gluon is a free parameter. Its 
precise value will depend on the UGD, on the assumed impact parameter distribution of valence charges in the 
nuclcon, on the dy/drj Jacobian and so on. 

In Fig. [5] we compare to ALICE and CMS data from refs. (43l - |46j . For UGD 111 we have adjusted the normal- 
ization relative to that determined from Pb+Pb collisions by a factor of 1.24; the energy and rapidity dependence 
of dN c h/dr] is consistent with the data. On the other hand, UGD 115 does not require any particular adjustment 
of normalization but appears to predict a slightly too steep energy dependence of the multiplicity and a narrower 
rapidity distribution (given our specific dy/drj Jacobian). 
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FIG. 4: Transverse momentum distribution of charged particles in the central region of p+p collisions at y/s = 7 TeV. The 
scale in the FF is taken to be Q = p±/2 (left) or Q = 2p± (right), respectively. 
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FIG. 5: Charged particle multiplicity as a function of pseudorapidity for p+p collisions at 
respectively, for two different UGD sets (glll9 and gllOl). ALICE and CMS data from refs. 



Js = 2.36 TeV and 7 TeV, 
43 4f| . 



V. MULTIPLICITY IN P+PB COLLISIONS 



In Fig. [S] we present our predictions for pj_-intcgratcd multiplicities in minimum bias p+Pb collisions. According 
to our findings from above, for UGD set glll9 we use the exact same normalization as for Pb+Pb and p+p 
collisions. On the other hand, we had found that the UGD set gllOl requires a correction of +24% to reproduce 
the measured multiplicity in p+p collisions at 2360 and 7000 GeV. Accordingly, for this UGD set and min. bias 
p+Pb collisions we have increased by hand the normalization by +10% (relative to Pb+Pb). A prediction for the 
charged multiplicity for the UGD with MV-model i.e. (7 = 1) can be found in ref . 1281 . We mention that the 
current results arc similar to other predictions based on the idea of gluon saturation p7H49| . This illustrates that 
indeed the energy and system size dependence of the multiplicity is governed mainly by the dependence of the 
saturation scale Q s on the target thickness and on x. 



VI. SINGLE INCLUSIVE SPECTRA IN P+PB COLLISIONS 



In this section we present single- inclusive charged hadron transverse momentum distributions for p+Pb collisions 
at y/ s = 5 TeV. For pseudo-rapidities near the central region, \rj\ < 2, the relevant light-cone momentum fractions 
in both projectile and target are comparable and small and so we use k± factorization, section IIII Al We fix the 
if-factor to the value extracted from p+p collisions as described in the previous section. Our default fragmentation 
function is KKP-LO evaluated at the scale Q 2 = kf. 
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FIG. 6: Predicted pseudo-rapidity distribution of charged particles in minimum-bias p+Pb collisions at y/s = 5 TeV for 
two different UGDs. 



We shall also show the "nuclear modification factor" i? p +pb defined as 



R P +Pb(p±) 



-Pb 



jdr] d 2 p ± 



dNP+ p /dr,d 2 P A 



(21) 



where (N co u) denotes the mean number of binary nucleon-nucleon collisions in a given centrality class; it is obtained 
from a standard MC Glauber model using an inelastic cross section of cri n (\/ s = 5 TeV) = 67 mb. For minimum 
bias collisions this leads to (N co u) « 7. 

To facilitate interpretation of the numerical results we shall not restrict to the AAMQS-likc UGDs with 7 > 1 
initial condition but also show some curves obtained with the UGD with 7=1 MV-model initial condition. We 
stress that although this UGD does not provide a good description of neither DIS data on protons nor of semi-hard 
Pt spectra in p+p collisions, it has not been directly tested against nuclear data yet and, therefore, remains a viable 
candidate for the initial condition for the evolution of nuclear wave functions. 
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FIG. 7: Average of the gluon densities of a single nucleon and a nine nucleon target divided by the gluon density of a five 
nucleon target, for Y = (MV i.e.) and Y — 3 (MV + rcBK), respectively. 

We first illustrate the effects of quantum evolution and of fluctuations in the thickness of the target. Fig. [7] 
compares the average UGD of a 1-nucleon and 9- nucleon target to that of a 5- nucleon target. The MV initial 
condition at Y = shows a strong suppression of this ratio at low intrinsic transverse momentum followed by a 
"Cronin-like" peak and an asymptotic approach to 1 from above (the leading higher-twist correction in the MV 
model is positive, ~ J rQ^/k^; appendix B in rcf. [50||). Neglecting fluctuations in the number of target nucleons 
could clearly distort the resulting R p a ratio significantly. Actually, this could be origin of the slight differences 
between our predictions for R p pb and those presented in [47] (shown in Fig.[TT]bclow) with similar dynamical input 
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but where the nuclear geometry is treated in a mean field approach, hence neglecting fluctuations. Rcsummation 
of small- a; quantum fluctuations removes the enhancement at intermediate kr- 




FIG. 8: Comparison of the rcBK-MC results obtained with only the elastic term of the hybrid formalism, Eq. (|19[) . to 
the RHIC forward data on single inclusive charged hadron (BRAHMS data [U) and neutral pion yields in p+p (left) and 
d+Au collisions (right). Solid lines correspond to the 7 = 1.119 i.e., dashed-dotted to also 7 = 1.119 i.c but using the 
prescription in Eq. (fTQ|) for the initial saturation scale. Dotted lines correspond to MV i.c. 

In Fig. E] we compare our results for single inclusive charged hadron (BRAHMS data [5l[ ) and neutral pion 
(STAR data (HI) distributions measured in p+p and d+Au collisions at RHIC. In this figure we include only the 
elastic component of the hybrid formalism. In what follows we adopt the DSS-NLO fragmentation functions as the 
default ones for all the calculations performed within the hybrid formalism. Our results show a good agreement 
with data. However, the figure also illustrates that RHIC forward data does not constrain well the initial conditions 
for the evolution of nuclear wave functions: both the UGD MV and glll9 sets (using either the natural, Eq. ©, 
or the modified, Eq. (|10[) . ansatz for the initial saturation scale at every point in the transverse plane) yield a 
comparably good description of data. This is due to the fact that transverse momentum distributions in the 
forward region do not probe the kx ^> Q s tails of the UGDs. 

Similar to previous phenomenological works, we found that no K-factors are needed to describe data at rapidities 
77 = 2.2 and 3. However, STAR data at more forward rapidities can only be well described if a K-factor rj 0.4 
is introduced. This may be an indication that large- a; phenomena non included in the CGC may be relevant in 
the region close to the kinematic limit of phase space. Note, however, that the value of the JiT-factor depends 
significantly both on the UGD and on the FF. 

In Fig. [9] we show the comparison to the same RHIC forward data, now also including the inelastic term in 
the hybrid formalism. We explore both fixed a s = 0.1 as well as one-loop running coupling at the scale Q. We 
observe that the effect of this additional term can be very large, especially at large transverse momentum. We 
note that, despite the fact that the coupling decreases with increasing transverse momentum, the running coupling 
prescription causes a larger effect than the fixed coupling one. 

We observe that the inelastic term exhibits a harder p^-dependence than the elastic contribution, and at some 
transverse momentum it overwhelms the elastic contribution. The crossing point depends on the particular choice 
of UGD. The effects from the inelastic corrections are stronger for the steeper glll9 initial conditions than for 
the MV ones over the entire range of transverse momentum shown in Fig. |9] Also, the importance of the inelastic 
term depends on the collision system or, cquivalently, on the target saturation scale: it is stronger for p+p than 
for d+Au collisions. For p+p collisions in particular it appears that the present formalism does not provide a 
stable result as the inelastic correction overwhelms the leading elastic contribution already at moderate values 
of transverse momentum. This not a completely unexpected result since, parametrically, the inelastic term is 
proportional to ln(p t /Q st ), with Q st the target saturation scale, while the elastic term scales as hi(p t /AQCD) (see 
discussion in (3lj]). Given the importance and magnitude of the inelastic term, our findings call for a complete 
phenomenological analysis of the full NLO corrections. 

We now proceed to p+Pb collisions at LHC energy, y/s = 5 TcV. In Figs. [TU1 and [TT1 we show our results for the 
single inclusive charged hadrons yields in p+p and minimum bias p+Pb collisions and the nuclear modification 
factor i? p+ pb for minimum bias collisions respectively. We compare also to i? p +pb from collinear factorization 
using EPS09 nPDFs [HD, [54| as well as to results from the "IP-sat" model and from an independent rcBK imple- 
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FIG. 9: Same as Fig. [8] but including the inelastic term in the hybrid formalism. Solid lines are the same as in Fig. [8] 
Dotted and dashed lines correspond to a s — 0.1 and a 3 = a s (Q = pt) in Eq. (|20[) . respectively. 



mentation 9 (47j . 

Before discussing the results let us first explain the meaning of the rcBK-MC bands shown in Figs. [TTIfTSl 
They comprise the results for i? p pb calculated according to Eq. (|2"Tj) using the three UGD sets (glll9, gllOl and 
MV), the three kind of fragmentation functions (KKP-LO, DSS-LO and DSS-NLO) and the two possibilities to 
determine the initial saturation scale (natural, Eq. ((9]), or modified, Eq. ([TU|) ) considered throughout this work, 
always using the same configuration in the numerator -p+Pb-spectrum- and denominator -p+p-spectrum-. The 
upper limit of the bands correspond in all cases to i? p pb calculated with UGD set gl.119 together with the 
natural prescription for the initial saturation scale. The black solid line in the plots for r/ = and 2 in Fig. [TT] 
represents the upper limit of the band if only modified initial conditions are used (such distinction is not necessary 
for Figs. [T2l and fT3l since both cases are treated separately). For the results obtained within the fct-factorization 
formalism the upper limit correspond to KKP-LO fragmentation functions, while for the results obtained within 
the hybrid formalism (both for only elastic and elastic +inelastic curves) the upper limit of the bands corresponds 
to DSS-NLO fragmentation functions. In turn, the lower limits of the bands correspond in all cases to UGD set 
MV and DSS-NLO fragmentation functions. The results for all other possible configurations - i.e. other UGDs 
and fragmentation functions and choice of natural or modified initial conditions- fall within the plotted bands; 
individual curves are not shown for clarity of the presentation. 



9 To mention two differences to our work: ref. |47|| uses a different fragmentation function and does not treat fluctuations of the 
nucleon configurations in the target. The predictions are not far apart but the difference illustrates the sensitivity of R p a to such 
"details" . 
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FIG. 10: Predictions for the single inclusive charged hadrons yields in p+p and minimum-bias p+Pb collisions at 5 TeV 
collision energy at rapidities 0, 2, 4 and 6. The grey bands at y=0 and 2 correspond to the rcBK-MC results using kt- 
factorization, Eq. (|13[1 . Dotted lines correspond to MV i.e.. In turn, the yellow bands at n — 2, 4 and 6 have been obtained 
using the LO hybrid formalism including only the elastic term, Eq. (|19|l . The overall normalization has been adjusted by a 
factor 50, 1, 0.02 and 4 • 10~ 4 for rapidities 0, 2, 4 and 6 respectively to improve the visibility. 



We observe a suppression of i? p +pb = 0.6 ±0.1 at midrapidity and pt = 1 GeV. Over the entire range of pt shown 
in the figure, and for both UGDs, i? p +pb decreases with increasing rapidity (towards the proton fragmentation 
region). This is consistent with expectations from non- linear evolution which suppresses additional gluon emissions 
in dense wave functions as compared to the dilute limit. However, we find that i? p+ pb increases with p t and reaches, 
or even exceeds, unity for hadron transverse momenta of several GeV. 

This is due to the following effects. First of all, the light-cone momentum fractions of gluons which contribute 
to hadron production at fixed rapidity increase proportional to p t . This drives us closer to the initial condition 
which in fact exhibits an "anti-shadowing" enhancement of the UGD tail at high intrinsic transverse momentum 10 . 
This higher-twist "anti-shadowing" component is further enhanced by Glauber fluctuations in the Pb target since 
<p(x,k±) is not linear in the target thickness while (N co \\) in the denominator of Eq. (|21[) is. Nevertheless, it is 
remarkable that small- x quantum evolution in the rcBK approximation predicts a disappearance at LHC energies 
(and r\ > 0) of the clear Cronin peak observed so far in all proton- nucleus collisions at lower energies 11 . Our 
results confirm earlier more schematic or qualitative predictions of this effect [5(| [5?} . Forthcoming LHC data will 
therefore provide a very important test for the evolution speed predicted by the running coupling BK equation. 

It is also remarkable that at intermediate rapidities rj sa 2 the hybrid formalism with only the elastic term gives 
-RpPb systematically below the one from fc t -factorization. However, the inclusion of the inelastic term Eq. ([20]) in 
the calculations within the hybrid formalism tends to systematically increase i? p pb as compared to only the elastic 
component, bringing it closer to the R p pb obtained from A: t -factorization. This effect is visible, in particular, at high 
transverse momentum, lifting R p pb closer to unity. This confirms the expectation that higher-order and energy 
conservation corrections should bring p t -spectra closer to the DGLAP limit at high transverse momenta, where 
gluon densities are smaller and non-linear corrections should be less relevant. However, as discussed previously, 
the precise quantitative effect of the inelastic term appears to depend quite strongly on the value of the strong 
coupling, preventing us to make precise quantitative predictions (only the results corresponding to a s = 0.1 are 
shown in Fig. (1111) : those corresponding to a s = a s (Q) produce a much larger effect on i? p pb). Also, it is somewhat 
surprising that the effect is strongest at the most forward rapidities r\ 6 and the highest transverse momentum 
(we cut off the curves corresponding to the inelastic term at pt ~ 8 in this plot, since they grow very fast above 
unity for larger p t ). However, this region (rj > 6, pt > 8 GeV) is close to the kinematic limit at LHC energies (see 
Fig. [5]). It could be that, similar to our discussion of forward RHIC data, the CGC formulation discussed in this 
work is not complete or inaccurate around the limit of phase space, hence the uncontrolled growth of i? p pb at the 



Note, however, that the UGD with 7 = 1 MV-model initial condition is very close to the EPS09 nPDF at high transverse momentum 
and y = 0. This indicates that for this UGD "anti-shadowing" is weak. 

1 Such a Cronin peak at 5 TeV is also predicted by some leading-twist shadowing models including parton intrinsic transverse 
momenta [55l |. 



16 




10 12 

pt (GeV/c) 



Rppb(r|=4) 

cme= 5 TeV 



rcBK-MC, min bias 
| rcBK-MC, Npart>10 
"j rcBK-MC, LO+inelastic term a=0.1 
] EPS09 nPDF 



1.5 




10 12 

pt (GeV/c) 



0.5 0.5 



Rppb(r|=6) 

cme= 5 TeV 



rcBK-MC, min bias 
' "; rcBK-MC, Npart>10 

| rcBK-MC, LO+inelastic term a=0.1 
~]EPS09nPDF 




1.5 



2 4 6 8 10 12 2 4 6 8 10 12 

pt (GeV/c) pt (GeV/c) 

FIG. 11: The nuclear modification factor i? p +pb for single inclusive charged hadrons in minimum- bias p+Pb collisions at 
5 TeV collision energy at rapidities 0, 2, 4 and 6. The grey bands at y=0 and 2 correspond to the rcBK-MC results using 
kt-factorization, Eq. I)13[l. In turn, the yellow bands at r\ — 2, 4 and 6 have been obtained using the LO hybrid formalism, 
Eq. (|19p . in minimum bias collisions. The blue bands between the dotted lines also correspond to LO hybrid results for 
collisions with a centrality cut iVpart > 10. Finally the dashed dotted curves at r\ = 2, 4 and 6 correspond to minimum bias 
collisions calculated within the hybrid formalism inch the inelastic term from Eq. (|20[) with ot s = 0.1. 



most forward rapidities. 

In Fig. [T2] wc show i? p +pb for two different centrality classes selected according to the number of participant 
nucleons 12 . At pt = 1 GeV we observe the expected pattern of stronger suppression (smaller i? p +pb) for more 
central collisions. In the -/Vp art > 10 centrality class suppression now persists up to p t — 2 — 3 GeV. 

For the UGD with 7 = 1 MV-model initial condition (lower end of the bands in Fig.[T"2"|) one observes, generically, 
the expected pattern: i) at y = there is suppression at low pt while i? p +pb — > 1 with increasing pt as the rapidity 
evolution window shrinks; ii) there is slightly stronger suppression at low p t for N palt > 10 central collisions while 
the centrality cut has very little effect at high p t ; iii) the suppression increases with rapidity and i? p +pb < 1 for 
all p t < 10 GeV at y = 2. 

The behavior of i? p+ pb with AAMQS UGDs (7 = 1.119 initial condition, upper end of the bands in Fig. [T2"|) in 
central collisions is more intricate. At p t = 1 GeV we still find the expected decrease of i? p +pb both with centrality 
and rapidity. However, for p t > 4 GeV we find that i? p +pb is very similar at y = and y = 2. This UGD exhibits 
rather non-linear (in the valence charge density) anti-shadowing at high intrinsic k t and so particle production at 
high p t in p+Pb collisions is dominated by fluctuations corresponding to a high valence charge density in the Pb 
target (high -/V part ). This can be seen from the fact that at y = 2 and high p t there is little difference between the 
minimum bias and iV part > 10 centrality classes. 



In p+A collisions it is not straightforward experimentally to perform centrality selection via impact parameter cuts. Also, because 
of large fluctuations impact parameter bins correspond to rather broad distributions of N pSjT t ■ 
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FIG. 12: Same as Fig. [TT] for two different centrality classes. Here, for all curves the initial saturation momentum squared 
in the Pb target is taken to be proportional to the density of nucleons per unit transverse area in a given event, Qs(xo; b) ~ 
Ta (b) , according to the natural prescription given by Eq. Q . Results in this plot have been calculated in the fc t -factorization 
formalism. 
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FIG. 13: Same as Fig.[l2]but here the initial saturation momentum squared in the Pb target is taken as Qg(xo; b) ~ T 1 ^ 1 
according to the modified prescription given by Eq. (|10[1 . Results in this plot have been calculated in the fa-factorization 
formalism. 



In Fig. [T3J finally, we show the nuclear modification factor for the UGDs with initial saturation momentum 
squared Q 2 s {xQ]h) ~ T^ 7 (b) which restores the Af(r;xo) ~ r 2 behavior of the dipole scattering amplitude in the 
perturbative limit, rQ s (xo) <C 1. This prescription for the initial saturation momentum reduces anti-shadowing 
at high intrinsic momenta. This leads to slightly lower values of i? p -|_pb at high p t than the AAMQS UGDs 
corresponding to Q 2 (xq\ b) ~ lA(b). 



VII. MULTIPLICITY AND TRANSVERSE ENERGY IN PB+PB COLLISIONS 

In this section we present the centrality dependence of the multiplicity and transverse energy in heavy-ion 
collisions at LHC energy 13 . This serves mainly as a rough check for the dependence of the saturation momentum 
on the thickness of a nucleus. We shall find that the present framework leads to a rather good description of the 
data. Nevertheless, we recall that we do not account for final-state effects such as entropy production; also, that 
our estimates rely on a crude "hadronization model" as well as on A: t -factorization which is not expected to be 



Very similar results were previously presented in l28l (unpublished) before the corresponding data was available 
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very accurate for p t integrated multiplicities in A+A collisions. That said, it is clearly justified to establish that 
the model is not in gross contradiction to basic features of heavy-ion data. 
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FIG. 14: Left: Centrality dependence of the charged particle multiplicity at midrapidity, 77 = 0; Pb+Pb collisions at 
2760 GeV. We compare our calculation for two UGDs to data by the ALICE collaboration. Right: Centrality dependence 
of the transverse energy at r\ — 0. 



We shall focus on the centrality dependence of the charged particle multiplicity at central rapidity, 77 = 0, which 
we determine along the lines described in section llll Al The result is shown in Fig.[14]for the UGDs with MV- model 
(7 = 1) and AAMQS (7 = 1.101 and Q 2 s (x ) ~ T A or Q 2 s (x ) ~ T 1 / 1 ) initial conditions. We use K = 1.43 for the 
former and K = 2.0, 2.3 for the latter. The number of final hadrons per gluon n g = 5 in all three cases. 

All UGDs give a rather similar centrality dependence of the multiplicity which is in good agreement with ALICE 
data [HI, On the other hand, they differ somewhat in their prediction for the transverse energy. This is of 
course due to the fact that the 7 = 1.1 initial condition suppresses the high-fc t tail of the UGD 14 . With the 
if-factors mentioned above these UGDs match the measured E t in peripheral collisions. This is a sensible result 
since one does not expect large final-state effects in very peripheral collisions. For the most central collisions the 
energy deposited initially at central rapi dity is about 0.5% of the energy of the beams and exceeds the preliminary 
measurements by ALICE [(30[ and CMS [Hf by roughly 50%. This leaves room for — p AV work due to longitudinal 
hydro expansion [62Tl63 |. 



VIII. CONCLUSIONS 



The upcoming p+Pb run at the LHC provides a new and unique opportunity to study the dynamics of very 
strong color fields in nuclei at high energies. The central question is whether QCD dynamically generates a 
semi-hard scale Q s which dominates particle production and to test its dependence on the valence charge density 
(resp. the nuclear thickness) and on energy. In this paper we have provided phenomenological predictions and 
expectations to the best of our ability to model the large- a; valence degrees of freedom paired with high-energy 
QCD evolution of the distribution of soft gluons in the rcBK approximation. As we have discussed at length 
throughout the paper, the CGC formalism at its present degree of accuracy is affected by sizable uncertainties. 
Some are due to the lack of quantitative control of higher order corrections to the formalism, while others are 
related to the lack of data constraining the non-perturbative parameters of the theory such as the initial conditions 
for the evolution or the impact parameter dependence of the UGDs. The forthcoming data will test whether the 
present formulation of the CGC effective theory is quantitatively accurate or not. 



One should keep in mind though that our estimate of the initial transverse energy carries a significant uncertainty of at least ±15% 
related to our choice of it-factor; it is not determined accurately by the multiplicity since the latter involves only the product of 
/^-factor and gluon — » hadron multiplication factor K g . 
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